This file evaluates seagrass gains and losses in coverage from 2018 to 2020. The goal is to characterize the distribution of depth estimates for areas where changes have occurred. A hypothesis is that seagrass may be lost at deeper depths as compared to where it was gained, possibly related to changes in water clarity that can cause seagrass to grow in more shallow areas where they are less light limited. Source content is here.
The analysis is in several steps:
First, the relevant libraries and datasets are imported. All analyses are conducted using the NAD83(HARN) / Florida West (ftUS) projection.
knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)
library(sf)
library(tidyverse)
library(tbeptools)
library(mapview)
library(units)
library(stars)
library(raster)
library(here)
data(dem)
data(sgseg)
data(sgdat2018)
data(sgdat2020)
# this is the projstring for NAD83(HARN) / Florida West (ftUS)
# https://spatialreference.org/ref/epsg/2882/
prj4 <- '+proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs'
# reproject layers
tbseg <- tbseg %>%
st_transform(crs = prj4)
sgmanagement <- sgmanagement %>%
st_transform(crs = prj4)
# reproject dem, convert to stars
demstr <- dem %>%
projectRaster(crs = prj4) %>%
st_as_stars()
# reproject seagrass layers
sgdat2018 <- sgdat2018 %>%
st_transform(crs = prj4)
stdat2020 <- sgdat2020 %>%
st_transform(crs = prj4)This code chunk does the following:
source as 2018 category, target as 2020 category)# 2018 seagrass data
a <- sgdat2018 %>%
left_join(fluccs, by = 'FLUCCSCODE') %>%
st_union(by_feature = TRUE) %>%
mutate(
Category = paste0(Category, ', 2018')
)
# 2020 seagrass data
b <- sgdat2020 %>%
st_union(by_feature = TRUE) %>%
left_join(fluccs,by = 'FLUCCSCODE') %>%
mutate(
Category = paste0(Category, ', 2020')
)
# so intersect doesnt complain about attributes
st_agr(a) = "constant"
st_agr(b) = "constant"
# union separate layers for faster intersect
aunion <- a %>%
st_union %>%
st_set_precision(1e5) %>%
st_make_valid() %>%
st_buffer(dist = 0)
bunion <- b %>%
st_union %>%
st_set_precision(1e5) %>%
st_make_valid() %>%
st_buffer(dist = 0)
# get full union
op1 <- st_difference(a, bunion)
op2 <- st_difference(b, aunion) %>%
rename(Category.1 = Category)
op3 <- st_intersection(a, b)
# combine and make multipolygon to polygon, do by id otherwise it doesn't work
unidat <- bind_rows(op1, op2, op3) %>%
mutate(
yr = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category))),
yr.1 = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category.1))),
Category = ifelse(is.na(Category), paste0('other, ', yr), as.character(Category)),
Category.1 = ifelse(is.na(Category.1), paste0('other, ', yr.1), as.character(Category.1)),
idval = 1:nrow(.)
) %>%
dplyr::select(idval, Category.1, Category) %>%
dplyr::select(idval, source = Category, target = Category.1) %>%
group_by(idval, source, target) %>%
nest() %>%
mutate(
data = purrr::map(data, st_cast, 'POLYGON')
) %>%
unnest('data') %>%
st_as_sf()The unidat dataset is subset to only gains and losses (i.e., polygons can remain the same between years or change between seagrass categories). Areas of each polygon are calculated in acres. Polygons less than 0.023 acres are also removed because these are smaller than the cell size (30 by 30 feet, or 0.023 acres) of the bathymetry layer. The file is saved because it takes a long time to calculate.
# select only lost/gained
chgdat <- unidat %>%
ungroup() %>%
filter(target == 'other, 2020' | source == 'other, 2018') %>%
mutate(
Acres = st_area(.),
Acres = set_units(Acres, 'acres'),
Acres = as.numeric(Acres),
var = case_when(
target == 'other, 2020' ~ 'lost',
source == 'other, 2018' ~ 'gained'
)
) %>%
dplyr::filter(Acres > 0.023) # remove slivers and those less than the pixel size (dem pixel size is about 30x30ft)
save(chgdat, file = here('data/chgdat.RData'))Reload the data for the rendered file. The results show polygon types in 2018 and 2020 in the source and target columns. A polygon type of ‘other’ in the source column means seagrass was gained in 2020 and a polygon type of ‘other’ in the target column means seagrass was lost in 2020.
load(file = here('data/chgdat.RData'))
chgdat## Simple feature collection with 4310 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 381530.6 ymin: 1150157 xmax: 531285.6 ymax: 1345319
## CRS: +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## # A tibble: 4,310 x 6
## idval source target geometry Acres var
## * <int> <chr> <chr> <POLYGON [US_survey_foot]> <dbl> <chr>
## 1 7 Patchy, 2018 other, 2020 ((411859.5 1250585, 411851.4 125~ 0.841 lost
## 2 8 Patchy, 2018 other, 2020 ((410664.3 1264177, 410671.8 126~ 0.237 lost
## 3 8 Patchy, 2018 other, 2020 ((411085.6 1262250, 411087.4 126~ 0.973 lost
## 4 8 Patchy, 2018 other, 2020 ((411220.3 1261246, 411225.4 126~ 0.239 lost
## 5 10 Patchy, 2018 other, 2020 ((391000.9 1273733, 390998.9 127~ 0.869 lost
## 6 15 Patchy, 2018 other, 2020 ((416303.4 1211170, 416301.8 121~ 15.0 lost
## 7 16 Patchy, 2018 other, 2020 ((416374.9 1210817, 416377.2 121~ 0.248 lost
## 8 17 Patchy, 2018 other, 2020 ((415834.4 1209104, 415845.7 120~ 0.902 lost
## 9 22 Patchy, 2018 other, 2020 ((385843.3 1274880, 385809.3 127~ 0.218 lost
## 10 25 Patchy, 2018 other, 2020 ((382545.8 1284993, 382537.6 128~ 0.109 lost
## # ... with 4,300 more rowsm <- mapview(dem, layer.name = 'Depth (m)') + mapview(chgdat, zcol = 'var', layer.name = 'Seagrass', col.regions = c("#F8766D", "#00BFC4"))
m